#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define N 250

float maior_valor(float *vetor, int ordem);


//-------------- funcaoes globais --------------------

//funcao sutracao de vetores
__global__ void add( float *a, float *b, float *c ) {

    int tid = threadIdx.x + blockIdx.x * blockDim.x;

        if(tid<N){
                c[tid] = a[tid] - b[tid];
        }
}


int main(){

    // --------- VARIAVEIS COM INFORMACOES DO ARQUIVO --------------------------
    FILE *fp;
    int ordem_matriz; //tamanho da matriz quadrada
    int fila_avaliada;
    double erro; //erro permitido
    int max_iteracoes; //numero maximo de iteracoes
    float **matriz_A;
    float *vetor_B;

    // -------- VARIAVEIS DO DEVICE
    float *dev_a, *dev_b, *dev_c;

    // --------- VARIAVEIS AUXILIARES ------------------------------------------
    int i,j,k;
    float **aux_matriz_A; // Diagonal de aux_matriz_A e nulo
    float *aux_vetor_B;
    float *aux_vetor_converge;
    int iteracao; // numero de iteracoes
    double parada_erro;
    float *x; // valores de x
    float *aux_x;
    float *diferenca_x;
    float aux,aux1;
    float aux2;
    float rowTest;

    // --------- LENDO INFORMACOES DO ARQUIVO ----------------------------------
    fp = fopen("matriz250.txt","r");
    if(!fp) printf("Erro ao abrir arquivo");

    rewind(fp);// posicao iniciao do arquivo
    fscanf(fp,"%d",&ordem_matriz);
    fscanf(fp,"%d",&fila_avaliada);
    fscanf(fp,"%lf",&erro);
    fscanf(fp,"%d",&max_iteracoes);

    fila_avaliada = fila_avaliada - 1;

    // aloca as linhas da matriz A
    matriz_A = (float **) malloc (ordem_matriz * sizeof(float *));

    // aloca as colunas da matriz A
    for ( i = 0; i < ordem_matriz; i++ )
        matriz_A[i] = (float*) malloc (ordem_matriz * sizeof(float));

    // aloca o vetor B
    vetor_B = (float *) malloc (ordem_matriz * sizeof(float *));


    // lendo matriz A
    for ( i = 0; i < ordem_matriz; i++ ){
        for ( j = 0; j < ordem_matriz; j++ ){

             fscanf(fp,"%f",&matriz_A[i][j]);

        }
    }

    // lendo matriz B
    for ( i = 0; i < ordem_matriz; i++ )
        fscanf(fp,"%f",&vetor_B[i]);

    fclose(fp);


    // allocate the memory on the GPU
    cudaMalloc( (void**)&dev_a, ordem_matriz * sizeof(float));
    cudaMalloc( (void**)&dev_b, ordem_matriz * sizeof(float));
    cudaMalloc( (void**)&dev_c, ordem_matriz * sizeof(float));

    // --------- VERIFICANDO SE EXISTE CONVERGENCIA ----------------------------

    // aloca as linhas da aux_matriz_A
    aux_matriz_A = (float **) malloc (ordem_matriz * sizeof(float *));

    // aloca as colunas da aux_matriz_A
    for ( i = 0; i < ordem_matriz; i++ )
        aux_matriz_A[i] = (float*) malloc (ordem_matriz * sizeof(float));

    // aloca o vetor B
    aux_vetor_B = (float *) malloc (ordem_matriz * sizeof(float *));

    aux_vetor_converge = (float *) malloc (ordem_matriz * sizeof(float *));

    //analisando convergencia nas linhas

    // aux_matriz_A recebendo seus valores
    for ( i = 0; i < ordem_matriz; i++ ){
        for ( j = 0; j < ordem_matriz; j++ ){

             if(i == j){ //se diagonal preenche com o valor 0

                  aux_matriz_A[i][j] = 0;
                  aux_vetor_B[i] = vetor_B[i]/matriz_A[i][i];
             }

             else aux_matriz_A[i][j] = matriz_A[i][j] /matriz_A[i][i]; // senao divide os valores da linha pela respectiva posicao i=j da mesma linha
        }
    }

    for ( i = 0; i < ordem_matriz; i++ ){
        aux_vetor_converge[i] = 0;
        for ( j = 0; j < ordem_matriz; j++ ){
            aux_vetor_converge[i] += aux_matriz_A[i][j];
        }
    }

    // analisando convergencia nas colunas
    if(maior_valor(aux_vetor_converge,ordem_matriz) > 1){

        printf("Nao converge...\n");
        exit(1);
    }

    // --------- CALCULA O NUMERO DE ITERACOES ATE A PARADA --------------------

    else{

    // aloca o vetor x
    x = (float *) malloc (ordem_matriz * sizeof(float *));

    // aloca o vetor aux_x
    aux_x = (float *) malloc (ordem_matriz * sizeof(float *));

    // aloca o vetor aux_valor
    diferenca_x = (float *) malloc (ordem_matriz * sizeof(float *));

    // inicializa os valores de vetor x
    for ( i = 0; i < ordem_matriz; i++ )
        x[i] = aux_vetor_B[i];

    iteracao = 0; //setando a variavel que armazena o numero de iteracoes
    parada_erro = 0.1;

    while(iteracao < max_iteracoes && parada_erro > erro){
    //while(iteracao < 3){

        for ( k = 0; k < ordem_matriz; k++ )
            aux_x[k] = x[k];

        for ( i = 0; i < ordem_matriz; i++ ){

            aux1 = 0;

            for ( j = 0; j < ordem_matriz; j++ ){

                  aux1 += aux_matriz_A[i][j] * aux_x[j];

            }

            aux = (-1) * aux1;
            x[i] = aux_vetor_B[i] + aux;

            // Utilizando CUDA na condicao de parada
            // copy the arrays 'a' and 'b' to the GPU
            cudaMemcpy( dev_a, aux_x, ordem_matriz *sizeof(float),cudaMemcpyHostToDevice);
            cudaMemcpy( dev_b, x, ordem_matriz *sizeof(float),cudaMemcpyHostToDevice);

            //diferenca_x[i] = aux_x[i] - x[i]; // diferenca de Xf - Xi
            add<<<ordem_matriz,1>>>(dev_a,dev_b,dev_c);

            // copy the array 'c' back from the GPU to the CPU
            cudaMemcpy(diferenca_x, dev_c, ordem_matriz *sizeof(float),cudaMemcpyDeviceToHost);


        }

        parada_erro = maior_valor(diferenca_x,ordem_matriz)/maior_valor(x,ordem_matriz);
        iteracao++;
    }

    rowTest = 0;
    for ( i = 0; i < ordem_matriz; i++ )
        rowTest += matriz_A[fila_avaliada][i] * x[i];

    printf("--------------------------------------------");
    printf("\n Iterations: %d",iteracao);
    printf("\n RowTest: %d => [%f] =?%f",fila_avaliada+1,rowTest,vetor_B[fila_avaliada]);
    printf("\n--------------------------------------------");

    printf("\n\nThe values of x[i] are, where i= 0,1,2,3... :\n\n");
    for ( i = 0; i < ordem_matriz; i++ )
        printf("x[%d] = %f  \n",i,x[i]);

    }
	
    return 0;

}

// FUNCAO QUE RETORNA O MAIOR VALOR DO VETOR
float maior_valor(float *vetor, int ordem){

      int i;
      float maior;
      float a,b;

      if(ordem == 1)
          return vetor[0];

      else if(ordem < 1){
           printf("\n\nErro ordem da matriz menor que 1...");
           exit(1);
      }

      else {
          maior = 0;
          for ( i = 0; i < ordem ; i++ ){

              a = fabs(vetor[i]);
              if(a>=maior)
                 maior = a;

          }
          return maior;
      }
}